

* Custom function
* ---------------
*
* my_fn is an arbitrary function with the following structure
*
* Required
* --------
*
* fn - existing variable to store the function value, f(x; a)
* fp - existing variable to store the first derivative, f'(x; a)
* x  - values where to compute the function
*
* Optional
* --------
*
* a - Additional parameters. This doesn't have to be called 'a', it can
*     be called anything (or you can omit it). All extraneous parameters
*     passed to my_nr are passed to my_fn.


capture program drop my_nr
program my_nr
    local kwargs MAXiter(int 1000) tol(real 1e-12) stol(real 1e-12) DISPevery(int -1)
    syntax newvarname, x0(varname) Function(str) [ `kwargs' *]

    tempvar x fn fp
    qui {
        gen double `x'  = `x0'
        gen double `fn' = .
        gen double `fp' = .
        `function', fn(`fn') fp(`fp') x(`x') `options'
    }

    local i  = 0
    qui sum `fn', meanonly
    local cmax  = max(abs(`r(min)'), abs(`r(max)'))
    local ctol  = abs(`cmax') > `tol'
    local cstep = 1
    while ( `ctol' & `cstep' & (`i++' < `maxiter') ) {
        qui {
            replace `x'  = `x' - `fn' / `fp'
            `function', fn(`fn') fp(`fp') x(`x') `options'
            sum `fn', meanonly
            local _cmax = max(abs(`r(min)'), abs(`r(max)'))
            local cstep = abs(`cmax' - `_cmax') > `stol'
            local cmax  = `_cmax'
            local ctol  = abs(`cmax') > `tol'
        }
        if ( mod(`i', `dispevery') == 0 ) {
            local mdisp `:disp %9.3g `cmax''
            disp "    iter `i', max `mdisp'"
        }
    }

    local --i
    local mdisp `:disp %9.3g `cmax''
    local tdisp `:disp %9.3g `tol''
    local sdisp `:disp %9.3g `stol''
    if ( `ctol' == 0 ) {
        disp "NR successfully converged:"
    }
    else if ( `cstep' == 0 ) {
        disp "No more improvements achievable:"
        disp "    step tolerance = `sdisp'"
    }
    else {
        disp "NR failed to converge:"
    }
        disp "    fun tolerance  = `tdisp'"
        disp "    iterations     = `i'"
        disp "    max |f(x, a)|  = `mdisp'"

    rename `x' `varlist'
end
